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I TK— VOLUME FORMULATION 


The integral form of the conservation laws for mass, momentum, and 
energy is given in equation (1). The solution vector 

W = [ p , pu, pv, pe] T , 


where p, (u,v), and E denote the density, Cartesian velocity 
components, and total internal energy, respectively. The quantities 
F and G are flux vectors. The variables x and y are Cartesian 

coordinates, and t represents time. In equation (1), 0 denotes the 
region of interest and dQ is the boundary curve. If the computa- 
tional region is partitioned with quadrilaterals and equation (1) 
is applied to each one, a system of ordinary differential equations 
is obtained. One of these equations is given in the figure. The 

quantity S is the area of the cell being considered. The vector W 
is evaluated at the cell center, and the indices i and j identify 
the cell. 


d_ 

at 


ff Wdxdy 

Q 


+ J (Fdy 

dQ 


Gdx) = 0 


( 1 ) 


Applying Eq. (1) to the quadrilateral ABCD 


d 

dt 




= 0 


where 

H AB’ H BC’ H CD' H DA '" f,uxes th rough sides of cell 


original page is 

OF POOR QUALITY 


FINITE-VOLUME MESH 


The fluxes ^AB * C ^ ^ CD y A 9 m the 

equations, account for momentum flow into 
convection and by molecular transport. If 


case t urn 

and of a* cell by 

the thin-layer Navier- 


Stokes equations are considered, viscous fluxes are included in only 


HfiC an d H da . Consider cell face BC. The viscous stresses at this 
face are determined by applying Green's theorem to cell A'B'C'D'. 


In a similar manner, the viscous stresses at cell face DA are 


obtained . 




NUMERICAL PROCEDURE 


In this ^figure^ the basic elements of the present procedure for the 
numerical solution" of the thin-layer Navier-St okes equations are 
given. A modified four-stage Runge-Kutta scheme is used to advance 
the solution in time. Artificial dissipation terms are added to the 
difference equations, and they are third order in smooth regions 
of the flow field. These terms are included for several reasons: 
(1) to enhance the coupling of the difference equations, (2) to 
control nonlinear instabilities, and (3) to eliminate oscillations 
at shock waves. Three techniques are employed to accelerate 
convergence to steady state. With local time stepping, the solution 
at any point in the domain is advanced at the maximum time step 
allowed by stability. This results in faster signal propagation and 
thus faster convergence. Enthalpy damping involves, in general, 
adding damping terms to the flow equations that are proportional to 
the difference between the local total enthalpy and the free-stream 
value. These terms vanish at steady state. Finally, residual 
smoothing is an implicit procedure to extend the local stability 
range . 


• Modified four-stage Runge-Kutta scheme 

• Artificial dissipation terms added (blending of 2nd and 4th 
differences) 

t Acceleration tech niques 

(1) Local time stepping 

(2) Enthalpy damping 

(3) Residual smoothing--extend stability range 


NUMERICAL PROCEDURE 


The computer code for the scheme is vectorized for efficient 
processing. The algebraic method of LeBalleur (ref. 1) is used to 
construct C-type meshes for airfoil calculations. The algebraic 
eddy viscosity model of Baldwin and Lomax (ref. 1) is employed for 
turbulence closure. 


• Computing efficiency enhancement -- vector ization of scheme 


• Algebraic mesh generation -- C-type for airfoils 


• Algebraic turbulence model -- Ba Idw in-Lomax 


i 
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FOUR-STAGE TIME-STEPPING SCHEME 


The four-stage scheme presented in the figure is used to integrate 
the flow equations in time. The quantities R, L, and D are spatial 
discretization operators. The operator L is associated with convec- 
tion, and the operator D is associated with the physical viscous 
terms and artificial dissipation terms. The terms of D operating on 
W are evaluated at the first stage and frozen for the remaining 
stages. For the given values of (i = 1,2,3), this scheme is 
second-order accurate in time for nonlinear problems. The advantage 
of the modified scheme over the standard Runge-Kutta scheme is that 
it requires less array storage in computer processing (an important 
consideration for three-dimensional problems). Central difference 
operators are employed for the spatial derivatives in the flow 
equations. The scheme is second-order accurate in space for 
sufficiently smooth grids. 


W (1) = W n - djAt R W n 
W (2) = W n - a 2 AtRW (1) 
w (3) = w n - a 3 AtRW <2) 

W (4) = W n + 1 = W n - AtRW (3> 
where on the (q + 1) st stage 

RW (q) « j (LW (q) - DW n ) 
and 

Oj = 1/4, a 2 = 1/3, a 3 = 1/2 
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ENTHALPY DAMPING 


In the present work, a damping term is added to the energy equation 
only. Note that in the absence of viscous terms, the total enthalpy 
( H) is constant throughout the flow field if there is a uniform free 
stream. In the case of viscous flows, constant H can be a solution 
of the energy equation. Moreover, if the dominant viscous terms 
(boundary-layer type approximations) are retained in the flow 
equations, the equation at the bottom of the figure is obtained. 
Then, if the Prandtl number is unity and the solid boundaries are 
adiabatic, constant H is a solution of the energy equation. With 
similar assumptions, the total enthalpy can be shown to be constant 
for turbulent flows. As indicated previously, enthalpy damping can 
accelerate convergence. It is especially helpful during the initial 
time steps of a calculation. At the present time, we are investi- 
gating other forms of damping so as to not preclude the capability 
to solve problems where heat-transfer effects are important. 


Consider the energy equation 


-|y (pE) + -|^ ( puH) + (pvH) + gp (_H - Hqq )" Viscous terms (1) 

damping term 


Retaining the dominant viscous terms, we obtain for laminar flow 


LHS (1) 


a_ 
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RESIDUAL SMOOTHING 


In this figure, some details of implicit residual smoothing are 
examined. A two-step scheme is considered. The quantity L is the 
spatial discretization operator. If residual smoothing is applied 
to the second stage, the factored form given in the figure is used. 
In this equation, the product symbol is taken over the number of 
spatial dimensions (£), 3 is the smoothing parameter, and <5 xx is the 

standard central difference operator. The quantities ^u n+ * - u n 
and |u^) - u n j are t ^ e updated and explicit residuals, respective- 
ly. If 6 is sufficiently large, the scheme is unconditionally 
stable. However, as revealed through linear stability analyis , the 
fastest convergence to steady state is not realized with a very 
large time step. 


Consider the two-step scheme 


(1) n n 

u * u - aAtLu 


u (2) = u n - AtLu (1) 


Residual smoothing for the second stage is given by 


TT 

i 





u 


n 


where 6 is central difference operator. 
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COMPUTATIONAL DOMAIN FOR AIRFOIL CALCULATIONS 


This figure shows a sketch of the computational domain for airfoil 
calculations. The outer boundary (ABC) is a C-type curve, and it 
reflects the type of mesh used. 


Outer Boundary 



Downstream Boundary 



BOUNDARY CONDITIONS 


During a calculation, a test is performed to determine if there is 
inflow or outflow at a given location along the outer boundary 
curve. If there is inflow, the total enthalpy H, entropy S, and 
tangential velocity component U are specified. If there is outflow, 
the Riemann invariant is specified. The quantities V, a, and y are 
the normal velocity component, speed of sound, and specific heat 
ratio, respectively. 


• SURFACE BOUNDARY 

(1) No slip 

(2) Adiabatic condition 

• OUTER BOUNDARY 

(1) I n f I ow --specify H, S, U 

(2) 0 u t f I ow --specify Riemann invariant 


• DOWNSTREAM BOUNDARY-- SPECIFY P 
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TURBULENT FLAT-PLATE FLOW 


The Runge-Kutta time-stepping scheme was used to compute turbulent 
flow (M ot = 0.5, Re w = 10 ) over a semi-infinite flat plate. In the 
figure on the left, the variation of the skin friction (based on 

free-stream conditions) with local Reynolds number is presented. 
The computed result is compared with an incompressible correlation. 
There is fairly good agreement except in the leading-edge region. 
Note that since the thin-layer Navier-S t okes equations were solved 
in the calculation, no attempt was made to resolve the leading-edge 

region of the plate. In the figure on the right, velocity profiles 

at X/L = 0.875 are compared. The present result agrees very well 

with that obtained with the implicit code of Steger (uses Beam and 
Warming scheme, ref, 1). 


(M = 0. 5, 
Correlation, M - 0 

OG 

C f = 0.0592 Re'- 2 
• x 

O NAVIER-STOKES 



Re * 10 6 ) 
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PARTIAL VIEW OF AIRFOIL GRID 


This figure shows a partial view of a representative C-type grid 
used in airfoil calculations. There is grid clustering in the 
leading- and t railing-edge regions of the airfoil. In addition, 
there is fine mesh spacing adjacent to the airfoil surface in order 
to resolve the boundary layer. 
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LAMINAR FLOW OVER AN NACA 0012 AIRFOIL 


This figure shows results for laminar flow over an NACA 0012 
airfoil. The free-stream Mach number and Reynolds number for this 
case are 0.5 and 5*10^, respectively. In the figure on the left, 
the variation of pressure coefficient with nond imens i ona 1 axial 
distance (X/C, where C is the airfoil chord) is presented. The 
viscous solution is compared with the inviscid solution. There 
are significant viscous effects in the leading- and trailing-edge 
regions. In the figure on the right, the computed sk in-f r i c t i on 
distribution is shown. Flow separation is indicated at the 

0.82 chord location. 


(M 


0.5, Re 


5xl0 3 , 


a = 0) 



x/c 
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LAMINAR AIRFOIL FLOW STREAMLINES 


This figure shows the streamlines on the upper trailing-edge surface 
of the airfoil. The indicated separation bubble is thin and 
elongated. Since the velocities inside the bubble are small, 
accuracy of the solution is of special concern. Therefore, two 

calculations were done. In the first case, the grid consisted of 
128 cells around the airfoil and 32 cells in the normal-like 
direction. In the second case, a 128x64 mesh was used. The cell 

spacing in the normal-like direction was half as large. The solu- 
tions were essentially the same except for the reverse flow 
region. This region was about 7 percent larger in longitudinal 
extent in the finer grid case. 



TURBULENT FLOW OVER AN NACA 0012 AIRFOIL 


Pressure and skin-friction distributions for turbulent flow over an 
NACA 0012 airfoil (M^ = 0.5, Re^ = 2.89x10^) are presented in this 
figure. There is good agreement between the computed and experi- 
mental pressure distributions. In the figure on the right, the 
skin-friction distribution calculated with the present Navier-S tokes 
code is compared with that determined with a boundary-layer code. 
The mesh spacing at the surface in the boundary-layer calculation 
was more than an order of magnitude smaller than that used in the 
Navier-St okes calculation. There is fairly good agreement between 
these solutions except in the immediate vicinity of the airfoil 
leading edge. The present solution probably does not have adequate 
resolution of the very thin turbulent boundary layer at the leading 
edge . 


(M - 0.5, Re = 2.89X10 6 , a - 0) 



x/c 
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SUBSONIC FLOW OVER AN NACA 0012 AIRFOIL 


In this figure, the predicted pressure distribution for subsonic 
viscous flow (M a, = 0.5 , Re*, = 2.93x10^) over an NACA 0012 airfoil 
at angle of attack (3.51°) is compared with that for inviscid flow 
and with experimental data. The angle of attack has been corrected 
for wind tunnel wall effects. There is good agreement between the 
present Navier-Stokes results and experiment. The influence of 
viscosity on this solution, in particular at the suction peak, is 
evident. 


(M - 0.5, Re = 2.93X10 6 , a - 3.51) 
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TRANSONIC FLOW OVER AN NACA 0012 AIRFOIL 


This figure presents pressure distributions for transonic flow 
( Moo = 0.756 , Rea, = 4.01*10^) over an NACA 0012 airfoil. The 
results from the viscous and inviscid calculations show good 
agreement with the experimental data. The shock in this case is 
very weak, and this is evident in the Mach number contours in the 
next figure. 


(M = 0.756, Re = 4. OlxlO 6 , a = 0) 



x/c 
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TRANSONIC FLOW OVER AN RAE 2822 AIRFOIL 


This figure shows the pressure distribution for supercritical 
flow over an RAE 2822 airfoil. The free-stream Mach number and 
Reynolds number are 0.730 and 6.5x10^, respectively. The wind 
tunnel corrected angle of attack is 2.79°. Transition of the 
boundary layer was fixed at the experimental value of 3-percent 
chord. In this calculation, a C-type grid with 260 cells around the 
airfoil and 96 cells in the normal-like direction was used. There 
is good agreement between the computed and experimental pressure 
distributions. The predicted lift and drag coefficients are also 
compared with those of the experiment. 


(M - 0.730, Re * 6.5xl0 6 , a - 2.79) 
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UPPER SURFACE SKIN FRICTION FOR RAE 2822 AIRFOIL 


This figure presents a comparison of the calculated and experimental 
skin-friction distributions. The sk in -f r i c t i on coefficient is based 
on boundary-layer edge conditions. In general, there is fairly 
good agreement with the experimental data. However, the computation 
shows a small shock-induced separation bubble and incipient 
trailing-edge separation, which are not indicated by the experi- 
mental data. The present computed result is generally consistent 
with those of other investigators. The predicted rise in skin 
friction at the airfoil trailing edge is probably due to inadequate 
representation by the basic Baldwin-Lomax algebraic eddy viscosity 
model of the transition from the wall boundary layer to the free shear 
layer . 


(M = 0.730, Re - 6.5xl0 6 , a - 2.79) 
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CONVERGENCE HISTORIES FOR LAMINAR AIRFOIL FLOW 


This figure displays variations of the common logarithm of the root 
mean square of the residual of the continity equation with the 
number of time steps. Curve A shows the residual decay when 
enthalpy damping is added to the Runge-Kutta time-stepping scheme. 
Curve B indicates the improvement in the convergence rate with the 
inclusion of implicit residual smoothing. Satisfactory convergence 
for engineering applications is achieved in 1,300 time steps, which 
corresponds to 2.5 minutes on the Vector Processing System (VPS) 32 
at Langley Research Center. 


(M 


0.5, Re 


5xl0 3 , a 


0, 128x32 GRID) 



0 500 1000 1500 2000 2500 
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A — enthalpy damping 


B — entha Ipy 
residual 


damping and 
smoothing 
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CONVERGENCE HISTORIES FOR TURBULENT AIRFOIL FLOW 


The figure on the left shows the convergence history when only 
local time stepping is employed to accelerate convergence. This 
convergence history is highly oscillatory. In the figure on the 
right, the convergence behavior for a computation using enthalpy 
damping and residual smoothing is much more monotonic. For this 
case, acceptable convergence is realized in about 1,250 time steps 
(4 minutes on the VPS-32). Note that the present computer code is 
not optimized for the VPS-32 system. With optimization the computer 
processing times can be reduced by a factor of 2 to 3. In conclu- 
sion, the final rate of reduction of the residual for the turbulent 
flow cases can be quite slow. Based on previous work with Euler 
equations, significant improvement in this rate appears to be 
possible with a multigrid scheme. At the present time, this is 

being investigated. 


(M ■ 0.5, Re * 2.89X10 6 , a ■= 0, 120x50 GRID) 


o 

D 

■o 

'(f) 

a> 

cr 


cn 

o 



130 


THIN-LAYER NAVIER- STORES (LAMINAR) SOLUTIONS FOR 
INFINITE RECTANGULAR WING 


In this figure, the velocity profiles computed from two- and three- 
dimensional thin-layer Navier-Stokes codes are compared at selected 
axial locations for an infinite, rectangular wing. The wing is 
constructed with constant chord NACA 0012 airfoil sections. 
Numerical solutions for this configuration were obtained at « 0.5 
and Re c = 5000. A 96x32 C-grid at each spanwise station was 
employed, and a total of four spanwise stations were used for the 
three-dimensional calculations. As expected, the solutions at each 
spanwise station were identical; and, as indicated by the results 
shown in this figure, these are in excellent agreement with the 
solutions obtained from two-dimensional analysis. 
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THIN-LAYER NAVI ER-S TORES ( TURBULENT) SOLUTIONS FOR 
INFINITE RECTANGULAR WING 


The configuration selected for the results shown in this figure is 
identical to the one shown in the previous figure. However, the flow 
here is assumed to be fully turbulent and the Reynolds number is taken 
to be 2.89x10^. A 120x50 C-grid was employed in each cross section 
for these calculations. Again, the results from the two- and three- 
dimensional analyses are found to be in excellent agreement. 
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PARTIAL VIEW OF GRID EMPLOYED FOR WING CALCULATIONS 


A partial view of a typical grid employed for wing calculations is 
shown in this figure. A C-H type grid topology is employed here; 
i.e., sections of C-type grids are stacked at various spanwise 
stations. Strong grid-clustering is employed near the wing surface 
to resolve extremely thin turbulent boundary layers that develop at 
high Reynolds numbers. The grid is also clustered near leading- and 
trailing-edge regions of the airfoil. A 121 x 51 x 17 grid is selected 
for an initial set of calculations. 
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CONCLUDING REMARKS 


A f ini te-volume scheme for numerical integration of the Euler 
equations has been extended to allow solution of the thin-layer 
Navier-Stokes equations in two and three dimensions. The extended 
algorithm, which is based on a class of four-stage Runge-Kutta time- 
stepping schemes, has been made numerically efficient through 
the following convergence acceleration techniques: (1) local time 

stepping, (2) enthalpy damping, and (3) residual smoothing. Also, 
the high degree of vectorization possible with the algorithm has 
yielded an efficient program for vector processors. The scheme 
has been evaluated by solving laminar and turbulent flows. Numer- 
ical results have compared well with either theoretical or other 
numerical solutions and/or experimental data. 


• Algorithm for Navier-Stokes equations constructed from a 
class of explicit multistage time-stepping schemes 

• Flexibility in treating arbitrary geometries obtained with 
a finite-volume formulation 

• Numerical efficiency achieved by acceleration techniques 

• Computer processing enhanced through vectorization 

• Scheme evaluated by solving laminar and turbulent flows 
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